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ABSTRACT 

Jordan canonical forms are used extensively in the literature on 
control systems. However, very few methods are available to compute 
them numerically. Most numerical methods compute a set of basis vectors 
in terms of which the given matrix is diagonalized when such a change 
of basis is possible. Here, a simple and efficient method is suggested 
for computing the Jordan canonical form and the corresponding trans- 
formation matrix. The method is based on the definition of a generalized 
eigenvector, and a natural extension of Gauss elimination techniques. 
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INTRODUCTION 

It is well-known that any matrix may be brought into the Jordan canonical 
form by a similarity transformation [1]. There are several methods available 
to compute the eigenvectors of a matrix when the eigenvalues are distinct 
[2-3], Some of these could be used to compute the eigenvectors for matrices 
with multiple roots. In Varah’s method [4] multiple eigenvalues are handled 
by perturbing the multiple eigenvalue to produce distinct eigenvalues. Eberlin 
and Boothroyd [5] also compute eigenvectors for multiple eigenvalues. However, 
none of these methods generate the basis vectors necessary to transform the 
given matrix into it’s Jordan canonical form. Chen [6] has suggested a 
procedure for computing the Jordan canonical form. Here, a simple and 
efficient algorithm, based on the notion of a generalized eigenvector, and 
using Gauss elimination techniques is given to compute the Jordan form of 
an nxn matrix. 


BACKGROUND 


- 1 , 


Given the nxn matrix A, we want to find the matrix T such that T AT 


is a Jordan matrix J. Let > A ) be the eigenvalues of A with multi- 
plicity ( n l» n 2 s » respectively. The number of eigenvectors 

associated with the eigenvalue A^ is given by = n-Rank (A-A^I) . The 
Jordan matrix, J, has the form 


J diag [J n , J 12 , , J. j J 21 , 3 1V » J 2a , I 

1 • / • 


9 





: j , , j 0 

. ml m2 


( 1 ) 
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with 


J ik = 


A. 1 

1 


o 


o 


1 


i = 1,2, ,m 


k = 1 s 2 5 } 


( 2 ) 


Let g., be the dimension of the block J.. and define 
ik ik 


i-1 £ k 

l B « + l V 

£— 1 j=l j“l J 


with a Q = 0 


(3) 


Let the generalized eigenvectors and the eigenvector corresponding to J ik 


be t 


, t 


— o./i n +l’ — a • r> 

i(k-l) i(k-l) 


• j t 


— o. ,, .+g.. -1 —a 

i(k-l) ik 


and t respectively. The 


ik 


transformation matrix T is made up of the n columns (t. , t~ , t , 

l i -o l;L 

-^ 1]L +1 S ^ 12 , .... +1 , ^ .The similarity 

l(a--l) la. ma 

1 1 m 

transformation satisfies the relation 


i.e. 


A [t.^ , s 


AT = TJ 




> t_2 * 


*9 t ]J 


-tr 


(4) 


Then, the eigenvectors of A f satisfy the relation 


(A-X r I )t ( , - 0 


£ - o o . , ..., o 
rl’ r2 


rot 


(5) 


Given an eigenvector of A^ the corresponding generalized eigenvectors satisfy 
the recursive relationship 


<A-» r « t*-! ' 1* 


1 °rk’ °rk _1 °rk _8 rk +1 


lc 1 j 2 , » • i j ct|^« 


(6) 
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The solution of equations (5) and (6) yields the transformation matrix T. 


COMPUTATION OF THE EIGENVECTORS ; 

Let A = (A - A^I) . We can choose non-singular matrices P r and Q r such 

that P AQ = U , where, U has the form 
r r r * r 

U = 
r 

Here is an (n-a^)x(n-a^) upper triangular matrix with l”ul * 0 and 
A^ is an (n-a )xa r matrix. Given (A-A^I) , P^, and U r can be obtained 
by Gauss elimination with full pivoting [7]. The eigenvectors corresponding 
to the eigenvalue A^ are obtained by solving the equation 

U t - 0 (7) 

r—£ — 

using a back substitution scheme employing a independent selections of 
the last a components of t_ 0 . Full pivoting guarantees that this will result 
in a linearly independent solutions which become the independent eigen- 
vectors corresponding to A . Substitution of these eigenvectors in equation 
(6) yields the set of generalized eigenvectors. 



Algorithm: 


1. Find the eigenvalues of A. Label them A^, A 2 , •••* 

2. Solve the equation U^t^ = £ for all eigenvectors corresponding to A r 
using independent selection of undetermined constants. The solution involves 


undefined variables v s w , ... . Generate an independent set of eigenvectors 

r r 

for A by setting each undefined variable in turn equal to 1 while holding 
r 

all other variables equal to 0. Denote the eigenvectors by t , t , ..., 

rl r2 
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3. For each eigenvector t ,i=l, 2 s ...,a form P Q and solve 

ri ri 

U t = P Q t 

r— a .-1 rr —a . 
ri n 

for generalized eigenvector corresponding to eigenvector t with the 

ri 

undetermined constants taking values given to them while evaluating t^ 

ri 

4. Repeat step 3 by forming P Q t , and solve Ut = P Q t .. 

r r-o ri -l r-o ri ~2 r r-a^-1. 

5. Continue to generate generalized eigenvectors as in step 4 until the 

equation U t =P Q t . becomes inconsistent i.e. when a non-zero 

r-o r-1-q r r— a -i 
ri J rx J 

quantity appears on the right hand side corresponding to zero rows of U^. 
This gives the basis vectors corresponding to the eigenvalue A r . 

6. Repeat step 2 thru 5 for r = 1, 2, . . . , m. to obtain all the basis 
vectors and hence the matrix T. 

7. Obtain the Jordan canonical form from J = T ^AT. Note that J need not 
be calculated directly since the block structure of (1) is determined by the 
number of generalized eigenvectors that are generated for each eigenvector. 


Computational Discussion; The computation of the eigenvectors and the 

generalized eigenvectors depend on the accuracy with which the eigenvalues 

of A are computed. Francis’ [8] algorithm is suggested for computing the 

eigenvalues. When the eigenvalues are approximate the calculation of the 

eigenvector can be refined as suggested by Wilkinson [9]. 

The algorithm suggested in this paper results in a large reduction in 

the amount of computation necessary to obtain the Jordan canonical form. 

til 

The number of computations necessary for an n order system with m distinct 
eigenvalues is shown in Table 1. 
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NUMBER OF COMPUTATIONS 


P 1 (A-A.I)Q 


I i + l i = T (n 3 -n) 


Total elimination 
for m eigenvalues 


m(n -n)/3 


Back substitution 


n-1 2 

1 I 1-^ 

i=l Z 


Total for n 
back substitution 


3 2 

< ILzS_ 

- 2 


To construct a 
right hand side 

(P^x) 


n 

l i - 

i=l 


2 

n -n 


Total R.H.S. 


0 3 2 

n(n -n)/2 = y- 


Total 


mn mn . 3 2 _ ,m+l 3 N 

~3 y + n -n =0 (— y n ) 
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A similar analysis of Chen’s algorithm [6] shows that the number of 
computations are of the order O^^). Thus the algorithm suggested here 
results in at least a fivefold saving in the number of computations. The 
method does not require the evaluation of the rank of matrices of powers 
of (A-X^I) as in Chen’s method. 


Examples ; 

The algorithm is applied to find the eigenvectors and the Jordan 
canonical form of two different matrices. 

A. Fourth order matrix: 

6 -3 4 1 

4 2 4 0 

4-231 
4 2 3 1 

This matrix is taken from Eberlin and Boothroyd [5]. The eigenvalues 
of the matrix are 5.23606797749979 (double root) and 0.763932022500210 
(double root) . 

The eigenvector and the generalized eigenvector associated with the double 
root 5.23606797749979 are 


0.4270509831 ' 


0.5868810394 ' 


1.0000000000 


1.0000000000 


0.3819660113 

and 

0.4721359550 

respectively 

_ 1.1458980340 


1.0901699410 _ 
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For the double root 0.763932G2250021Q the corresponding vectors are given 
by 


’-0.3726779962 


’ 0.2197175016 ' 

0.1273220038 

and 

0.4182146692 

0.3333333333 


-0.3171224407 

1.0000000000 


1.0000000000 


Notice that the two eigenvectors and the two generalized eigenvectors are 
all independent unlike in [5] . The Jordan canonical form can be readily 
written as 


5.2360 1 i 

i 

0 5.2360 j 

o 

i 

i i 

0.7639 1 

O | 


i 

i 

0 0.7639 


The execution time was 1.57 secs with a WATFIV (Univ. of Waterloo - Fast 
Fortran) compiler. 

B. System matrix of Boeing Helicopter 

The following 8x8 matrix arises in the design of a helicopter 
stabilization system using Pole-placement theory [10] . 
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.021 

.025 

-29.64 

.6968 

.1879 

0 

-.0941 

0 

-.0903 

-.802 

-80.98 

-1.878 

.5524 

0 

-8.517 

0 

0 

0 

0 

1 

0 

0 

0 

0 

-.0058 

.0145 

1.4672 

-1.460 

.45 

0 

.068 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

-784 

-35 

0 

0 

0 

o 

0 

0 

A 

0 

0 

1 

0 

0 

0 

0 

0 

0 

-784 

-35 


The eigenvalues of the system computed by using Francis' method are 
0.50432908, -2.3585084, -0.19350035 + j 0.35283477 and -17.5 + j 21.857493 
(double root) . The eigenvectors corresponding to the distinct roots are 


’ 1.000C000C00 


0.2528902161 


-0.0949009676 + j 0.6460398691 

0.9167473189 


1.0000000000 


1.0000000000 + j 0.0000000000 

-0.0157197678 


G.G2GC347219 


-0.0074706563 + j 0.0035914411 

-0.0079269851 


-0.0472520599 


0.0026856551 + j 0.001S5Q1968 

0.0000000000 


0.0000000000 


O.OGOOOCOOGO + j 0.0000000000 

0.0000000000 


0.0000000000 


0.0000000000 + j 0.0000000000 

0. 0000000000 


o.oooooooooe 


0.0000000000 + j 0.0000000000 

0.0000000000 

9 

. O.OOOQOOGCOO 

9 

. 0.0000000000 + j 0.0000000000 _ 


respectively. Each of the double roots has two eigenvectors associated with 


it. These are 
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-0.0000183498 + j 0.0002379966' 


' 0.0000224177 + j 0.0001119734 

-0.0001564383 + j 0.0007421192 


0.0026667158 + j 0.0107258554 

0.0000193897 + j 0.0000084539 


0.0000031152 + j 0.0000011743 

-0.0001545381 + j 0.0005715717 


-0.0000288496 + j 0.0000886422 

-0.0223214285 + j 0.0278794553 

and 

0.0000000000 + j 0.0000000000 

1.0000000000 + j 0.0000000000 


0.0000000000 + j 0.0000000000 

0.0000000000 + j 0.0000000000 


-0.0223214285 + j 0.0278794553 

0.0000000000 + j 0.0000000000 


.. 1.00000000C0 + j 0.0000000000 


Since the multiple eigenvalues have as many eigenvectors as their multiplicity, 
the Jordan canonical form for this matrix is diagonal and is given by 

diag [.50432908, -2.3585084, -0.19350035 + j 0.35283477, 

-0.19350035 - j 0.35283477 -17.5 + j 21.857493 -17.5 + j 21.857493 

-17.5 - j 21.857493, -17.5 - j 21.857493] 

The execution time using a WATFIV compiler was 8.69 secs. 

Flowchart and Computer Program ; 

These are given in Appendix A and Appendix B, respectively. 

Conclusion ; 

A method has been outlined to find the basis vectors to transform 
a given nxn matrix to its Jordan canonical form. The method is simple 
and efficient. It does not require the evaluation of the rank of matrices of 
powers of (A-X_^I) as in Chen’s method [6]. There is at least a fivefold 
reduction in the number of computations. Two examples are given to 


illustrate the method. 
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APPENDIX A 


\ 

^ CANON J 


Read NxN Matrix I 

A and NE | 

distinct eigenvaluesj 
— x 


L=1 


Form (A-A £ I) 

f 


CALL UMAKER 


/* Lu Decomposition 
of (A-A £ I) 


Test UM for 
Rank and # of 
eigenvectors 


Call EIGVEC 


/Print eigenvector 
and j 

generalized / 

eigenvectors . / 


/* Eigenvectors & generalized 
eigenvectors corresponding 
to X £ */ 



STOP 




Main Program 







A-2 



/* Preset undetermined 
constants for iith eigenvector 
*/ 


Back Substitute 
to find eigenvector 


4 


YES 


Cntcol>n 

!no 


/Too Many 


I 


eigenvectors / 


STOP 


Fill in NTCOL of j 
Modal matrix and ! 
form the ne w R. H.S. i 

\ 

J Operate on j 

| by row operations j 

| NTC OL =NT C OL+1 
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YES 


r 

i 


V 














B-l 


1 

„ Z. 

3 


4 

5 

6 

7 

8 
0 

10 

11 

12 

13 

14 

15 

16 

17 

18 


IS 

20 


APPENDIX B 


C MAIN PROGRAM 

C 

I M°LIC IT PEAL *3(A-F,0-Y) 

COMMON AM.PMB, ,PM! ,IJMR,UMI ,TMR,TMI,PEIGP*PE1GI»PMER,PMEI»YR* Y.I , 

1EPSA,E?R, IR* 1C* I E IG»NTCOL *LL ,ME»M.Nfc IG.IRANK . I OPT 
01 M£M5 1 AM ( 12 » 12 )* PMS ( 12 * 12) « °M1( 12* 12 ) *U-‘.R ( 12, 12) , 

1UM!( 12, 12) ,TMR( 12, 121 ,TMI< 12,121 ,PEIGR<12) , PE I Gil 1 2 1 , PMER (12 , 12 1, 
2PMEI(12»12)»YR(12)»Y1(121»IR(12)*IC( 12) » I E I GI 15,16) 



THIS F0UT1NE IS OESIGNEO TO FIND ALL THE EIGENVECTORS ANC 
GENERALIZED EIGENVECTORS OF AN REAL MATRIX GIVEN THE SET 

OF DISTINCT EIGENVALUES. THE PPINTCUT INDICATES THE 
APPROPRIATE JO°OAN_. CANONICAL FCSM - 

VAR! AP-LES 

ALL VAR I AGUES FROM A TO H ARE DOUBLE PRECISION 

ALL VARIABLES FROM 0 TO Y ARE COMPLEX EI T F SEAL PART 

ENDING IN S AMD IMAGINARY PAST ENDING IN I 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


AM ORIGINAL MATRIX 

PE I G ,R ♦ I EIGENVALUES. __ - 

PM, R*I < A -LAMBDA* I) MATRIX 

UM.PtI DECOMPOSED MATRIX, U ABCVE DIAGONAL, M BELOW 
TM,R+I MATR I X OF GENERALI ZED E IGENVECTORS— MODAL MATRIX _ 

N 01 MENS ION CF AM 

NE NUMBER of eigenvalues - 

IC " CCL'JMN INTERCHANGE INDEX VECTOR 

IP ROW INTERCHANGE INDEX VECTOR 

IOPT 0 PTI 0 N< = 1 ) FOR INTRMEDIATE PRINTOUTS 

I END OPriONl = l) FOP ADDITIONAL PROBLEM TO FOLLOW 

4000 FORMAT (IHl) 

4001 FORMAT! 81 10 ) 

4002 FORMAT ( 4020 . 10 ) .. - 

4003 F C c MA T ( /////) 

4004 FORMAT!//) 

4005 FORMAT ( 2030 . 15 ) ... 

4010 FORMAT! //, 5 X, 'Mi TR IX DIMENSION = ', 13 , 

1 1 number of distinct eigenvalues = ' 13 , /) 

4011 FO c MAT ( 5 X, • A MATRIX',/) 

4012 FORMA 7 ( 5 X , ' 0 1 ST I NC T EIGENVALUES' » // » 5 X » 9 H? EZL PART , 9 X , 

114 F IM\G INARY PART,/) 

4013 F )F MA T( 5 X , ' R CW INTERCHANGE INDEX',/) ... 

4014 FORMAT ( 5 X, ' COLUMN INTERCHANGE INDEX',/) 

4015 FORMAT ( 5 X. • CECOMPOSED MATRIX',/) 

4016 FORMAT ( 5 X ,' MU w BfcP OF EIGENVECTORS CORRESPONDING TO EIG ', 12 , 

O' IS • , I 3 , / ) 

4017 FORMA 1 ' ( 5 X ,' BLOCK NUMBER ', 12 ,' FAS ', 12 , 

1 ' GENERALI ZED E I GENVEC TOP. S' . / , 5 X , • THE FIRST IS THE EIGENVECTOR’,/) 

4018 FORMAT ( 5 X, OHREAL PART, 9 X, 14 HIMAG INAP Y PART,/) 


C 

C _ INPUT A MATRIX AND EIGENVALUES 
r 

10 READ 4001 ,\ ,NE » IOPT ,! END 

PRINT 4C00 
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21 

22 

23 

29 

25 

7 A 


PRINT 9003 
PRINT 9010, M,NE 

READ 9002, C 1 AMI I ,J ) »J= 1,N 1 » 1 = 1,N 1 
PRINT 9011 

PRINT 900?, 1 I AMI I » J I , J*1,NI , I a l ,N) 

DC 50 1 = 1, N - 

to 

27 

. 

DO 50 J=1,N 

28 


TMR II , J ) = O.COO 

29 


TMI II , J > =0 .000 

30 


PMR 1 I , J »=AM( I , J) 

31 

50 

PM II 1, J 1 = 0 .000 

32 . .. 


PRIN T 6009 - 

22 


PRINT 9012 

39 


DC 60 1=1, NE 

35 


READ 9005,4,8 . .. 

36 


PRINT 9005, A, 0 

37 


PFIGRII ) = A 

3e 

60 

PEIGIin=B 


c 

c 

c 

39 

90 

91 . 

92 

93 


START COK.PUT IMG 

PR TNT 9000 
PRINT 9003 

NTC 0L = L 

DO 500 t.= l»NE 
LL=L 


C 

C 

C 

99 

95 

96 

97 

98 

99 

50 

51 

C 

S 2 
53 

C 

C 

C 

59 

55 

56 

57 

58 

59 

60 
61 
62 

C 

C 

c 

63 

69 

65 

66 
67 


FORM ( A- LAMBDA* I 1 

DO 100 I = 1 * N 

DO 100 J=1*M 

PM E p ( I , Jl = PMR ( I , J 1 

P«FI C I » J) = PMI ( I .Jl 

IF ( I- J >100,95, IOC 
95 pME^it.jjsP^EPItfJl-PEIGPIL) 

pm Fill , J)=PMEI < I ♦JI-PF.IG! C.L) 

100 CONTINUE 


CALL UMAKER ----- 

I F ( IOPT 1 107, 107, 105 

IF IOPT=l PRINTOUT DECOMPOSED MATRIX 


105 PRINT 
PP I NT 
PRINT 
PRINT 
PRINT 
PRINT 
PR !■* IT 
dp ;\'T 
PP I NT 


9003 
9 013 

9001, ( IRIII, 1=1, N> 


9019 

900 l, ( IC( I ),I = 1,N» 

9015 

9002,1 (UPS (I ,J) »I 


9009 

9002, ( IUMI lt,J) ,J=l,N),I=l,N) 


TEST UM FOR RANK AND NUMBER OF EIGENVECTORS 


107 CONTINUE 

CALL CABS(EPSA,UMR(N,N>,yMI(N,N)l . ... 

EPS A=E PSA* 100 .000 
IF(EPSA-l.0O- 12) 11C, 110, 112 
110 £ PF= l , OG- 1 2 - 
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68 CO TO 115 

69 112 EPR=EPSA 

70 115 CONTINUE 

71 NE IG= 1 

72 NM = f|-l 

,_13 DC 125 I=i»NM 

74 I A=fl— I 

75 CALL CABS! ATEST,UM*UA,IA),UMII IA, IA) ) 

76 IF(ATEST-E?C )120, 120*118 

77 118 GO TO 130 

78 120 NE IG=N£ IG*l 

75 125 CONTINUE . . 

80 130 CONTINUE 

81 lEIG(L.l) =NE IG 

8Z PRINT 4016»L*NE IG 

C 

63 CALL 6IGVEC 

: C ~ PRINTOUT RESULTS 

C 

. 84 .. . IC T = l : 

85 DC 135 KKK = 1 » L 

£6 IFIKKK 1)132,135,135 

87 132 KP = IE IGIKKK, l )♦ 1_ 

88 DC 133 KC = 2« KB 

85 133 1C T=ICT*IEIG{ KKK ,KC 1 

90 135 CONTINUE 

- 91 00 150 J=l,NE!G 

52 JJ=IEIG(L,J*1) 

93 PR I NT 4C04 

94 PRINT 4017, J,JJ 

55 PRINT *018 

56 DC 150 K = l,JJ 

57 PRINT 4004 

93 00 140 K K = 1 , N 

99 PRINT 4002, T M P(KK, ICT), TNI ( KK » IC T I 

100 140 CGNTINUE 

1C l ICT-ICT+1 

102 150 CONTINUE . .. 

C 

103 500 CONTINUE 

f 

C TERMINATION" 

C 

104 IF ( IE NO 1510,510, 10 

1C5 5 10 CONTINUE 

106 PRINT 4000 

107 STOP . ... 

1C8 END 


5t)as-ni)T piE U M AKER 

IMPLICIT REAL *8 ( A -H, C-Y 1 

COMMON AM , PMR , PM I , U V R . , IJMI , TMR , TMI , PE ICR , PE IG I , PMER , PME I ,YR, Y I , 

IE PS A, EPF , IM, IC, I E I G, NTCOL , LL , N E , N , NE I G , I RANK , l OPT 
01 MEN SIGN am (12, 12) ,PMR (12,121 , I! 12, 12 1,'JMR ( 12,121 , 

Him [ { 12, 12) ,7MR( 12, 12) , TMI ( 12,121 , PE IGF ( 12) , PE I G It 1 2 1 , PMER ( 12 , 12 1 , 
2PME M 12 , 12 1 , YRl 12 I , Y I ( 12 1 , IR( 12 1 , IC C 12) , IEIGI 15, 16) 

THIS SUBROUTINE CALCULATES THE LL' DECOMPOSITION CF 
P M £ , F. ♦ I WITH FULL PIVOTING 


109 

no 

in 

112 
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C 

c 

c 

. ___.c. 

c 

c 

113 . 

114 

115 

116 

117 

118 

C 

C 

c 


115 

120 

121 

C 

c 

c 

122 

123 

124 

125 

126 

127 

128 
129 

C 

C 


120 

131 

132 

133 

C 

c 

c 

124 

135 

136 

137 

138 

125 

C 

C 

c 

140 

141 

142 


INPUT VARIABLES 

PM£,R*l MATRIX TO ee DECOMPOSED 
N DIMENSION 


OUTPUT VARIABLES __ -- •• 

UM,R,I DECOMPOSED MATRIX, L UPPER TP I ANGLE INCLUDING 
DIAGONAL, multipliers felow diagonal 

ip _ pnv» interchange, index vectcr, 

IC COLUMN interchange index vector 


preset row and column interchanges 

DO 100 1*1, N.. 

IP 1 1 >= l 

ioo icm = t 

EP S = l • OD- 2 0_ , 

IP ANK = N 
NN=N-l 

BEGIN E L I M I NAT I ON P ROC E CUR E 

00 200 LS= l ,NN 

LSS=LS 

AMAG=0.0D0 

SEARCH FOR PIVOT 

DO 120 I =L S ,N , ' , - 

DO 120 J=L$,N 

CALL CAfiSI ACUM, PMER I IR ( I 1 , IC ( J 1 ) , PMEI ( IR 1 1 1 • I C ( J ) ) 1 

I F ( AD’JM- AMAG 1 120 , 120,105.. .... 

105 IS=I 
J$ = J 

AM AG=AOL M -- 

120 CONTINUE 

TEST FOR COMPLETION 


IFI AM AG- EPS 1125, 125,130 

125 I RANK =L S- 1 _ 

GO TO 300 
13C CONTINUE 

INTERCHANGING ROW AND COLUMN INDICES 


I T = !P ( LSI ... . ... 

IP(LS)=lRnS» 

IP ( 1S1 =IT 

IT = IC ( L S I - 

IC ( LS 1 = IC I JS 1 
ICIJSMIT 

ELIMINATE IC (L S » COLUMN AND SET UP UM,R*I LS ROW 

L SP = L S 1 _ . - - - 

DO 150 I=ISP,N . 

CALL CU!V(QP«C!,P M CMI*<(I)»IC(LS)),PMEI<IR(I)»IC(LS))» 
1PMER 1 1 R I LS ) , IC ( L S 1 1 , PMfc l IIP (LSI , 1C I LS 11 1 
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143 

144 

145 

146 

147 

148 

149 


UMS( I ,CS) =CR 
UMI<I.LS)=Of 

CALL^MUL ( CTR.CTI .CF.QI ,PMER» !R<LS ). IC 1 J I) . PME I ( IR1LS 1 , IC ( J 111 
PM£R< Ull l,IC( JI»=PM6fi(IRU).ICC J) ) OTR 

PHEI < JR ( 1 ) i JC.(.J) )=P y EIJ.|R<l. )* IC JJ UrQTL 

ISO CONTINUE 


C 

C . .PATCH UP. RANK TEST .._ 

C 

150 1F1LS-NN) 170, l£C, 16C 

151 __ 160 CALC C A eS < AC'JM , PM ER ( 1* C N> , IC ( N I 1 » .PME. I ( I R < N ) , I C I N ) ).) 

• 152 I F < ADUM-EPS ) 165 » 165 .170 

153 165 I«AMK=M-1 

154 170 CONTINUE 

155 200 CONTINUE 

156 GO TO 350 


C WINDUP PROCEDURE 

C 

157 300 00 310 I =1 f N 

158 00 310 J=t .N 

159 UPP(I.Jl=P M ER(IRm» ICIJll 

14C _ UMK I . J) *PRE! ( IP! ! ) » ICC J)_) 

161 ' ’ IF<I-J)302,31C,310. 

162 302 IF(I-LS)310,304,304 

163 304 UM? ( J » I ) =0. 000 

164 UM M J » I 1 = 0*000 

165 310 CONTINUE 

166 GCJ TO 400 ... ..... ... — 

167 350 DO 360 I»l,N 

168 00 36)0 J 3 1 » N 

165 u ms ii ,j)=PMFR<iR(n .icon 

170 UMI (I , J)=PME! IIRI I)» ICI J) ) 

171 360 CONTINUE 

172 400 COf! T l'IUE .. 

173 RETUPN 

174 EMO 


175 

176 

177 

178 


179 


C 

C 

c 

r 

r 

C 

C 

c 

c 

c 

c 

c 

c 


SUBROUTINE eigvec 
IMPLICIT NEAL *5Ta-F,0-Y) 

COMMON fcM , PMP , PMl ,UM' J »UM|»TMR»TMI,PEIGR,PEIGI»PMER,PMEr»YR»YI, 
1EPS A, EPP » I s , IC. IE IG.N7CDL .Ll.NE .N.NEIG.IRANK, 10PT 

DIMENSION AMI 12 ,12) .PMR 112. 12 ) ,°MI ( 12.12 J.UMRl 12. 12) , 

IUMM 12, 12 » ,7MR( 12,12) *TMI I 12,121 ,PriGRCl2) , PE IG I ( 12 ) , PMER <12 , 12 ) , 

2PME I <l2»12)»YRli2)»YHl2),iRIi2)»IC(l2)»IElG(l5»16) 

01 MEN SI ON RP ( 12) .RI <12) »SA( 12) ,SB(12) »IRA( 12) 

THIS SUBROUTINE TAKES THE DECOMPOSED MATRIX OF UMAKER WITH 
KNOWN CAN* ( f’i-NE 1 0) ANO CALCULATES AIL THE EIGENVECTORS AND 
GENERAL I ZEO EIGENVECTORS OF THE CURRENT EIGENVALUE <P£IG(d ) 

INPUT VARIABLES 

L)M, R ♦ I DECOMPOSED MATRIX 

N DIMENSION 

IR.JC ROW AND COLUMN INTERCHANGE INDICES 

NE1G HUM PE 5 0 F EIGENVECTORS _ ... .... 

NTCOL CURRENT COLUMN CF TM 


OUTPUT VA° I ABC ES 



c 

c 

c 

c 

c 

c 

c 

c 


TH,RH COLUMNS OF “HDAL MATRIX - ALSO EIGENVECTORS 
AND GENERALIZED EIGENVECTORS 

IEIG NUMBERS OF EIGENVECTORS AND GENERALIZED EIFENVECTOPS 

CORRESPONDING TO EACH EIGENVALUE 


BEGIN SEARCH FOR EIGENVECTORS 


1E0 NOK=N-NE IG 

181 00 200 11=1, MEIC 

182 NUP = l 

1E3 _ N I =N«-1- 1 1 _ _ 

C 

C PRESET UNCETERM INED CONSTANTS FOR II--TH EIGENVECTOR 

18 A DO 50 J = 1,N 

185 RR I J ) =0 .000 

lEt PI(J) = 0.000 _ 

187 YR<J)=0.000 

188 50 YI ( J) =0.000 

189 YR ( MI ) = t . 000 _ _ . . .. 


C 

C BACK SUBSTITUTE TO FIND EIGENVECTOR 

C _ _ 

190" 60 CONTINUE 

151 DO 75 J = 1 , NOK 

192 J J = NOK 4-1- J _ 

193 JK=J»NEIG-l 


199 

155 

196 

157 

158 

199 

200 
201 
202 
203 
2C4 


70 


DO 70 <=1.JK 
KK =N+1 -K 

CALL CMiJLIOTR.OTI ,iJMR( JJ,KK J ,UNI( J j ,KK ) , YR I KK » , VII KK ) » 
SAIKI = -CTR 

SB I K ) = -QT I 

CONTINUE 

CALL SIY( JK,SA,SMR) 

CALL 3'JM( JK.SB.SMII 

SKC=SMR+PR (JJ I 
SMI=SM1+RI( JJJ 


205 



YR ( J J J =CTR 



2C6 



YI( JJ)=OTI 



207 


75 

CONTINUE 

- - 

■ 


V 

C 


FINO ALL GENERALIZED 

EIGENVECTORS 


C 





208 



NGE = l 



209 


76 

TFCNTCUL* U» 75,79,77 



210 


77 

PRINT 605G 



211 


5050 

FOPMAT ( 5X , ' TCC PANY EIGENVECTORS 

FOUND' ,// ) 

212 



STOP 



213 


79 

DC 80 ! = 1 , N 



215 



TP?( IC« I) ,NTCOL» =YR(I) 



215 

c 

80 

TM I ( ICI n,NTCOL) = YHI) 




c 


OPERATE CN RIGHT HANC 

SIDE BY 

PCW CPS 

216 



DO 90 1=1, N 



217 



ioAm=i 



218 



RR < 1 )=TMR( I.NTCOL) 



219 


90 

RI U ) = ^HI (I.nTCOL » 
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220 

221 

222 

223 

224 

2 2.5 

226 

227 

228 

229 

23C 

211 

232 

232 

234 . . 

235 

236 

2 37 . .. 

238 

239 

240 

241 

242 

243 

244 

C 

C 

c 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 .... 

256 

257 
256 

259 

260 


261 
262 
263 
26 4 


NTCOL =NTCOL* 1 
NP = N- 1 

00 120 1*1. MM 
IF ( IP.AI t 1-1-1 1)194, ICC, 94 
94 00 93 ! S= I »N 

I F < I - A < I S » r I R « l ) > 9 8 , 96 , 98 

96 isr= is 

98 CONTI HUE 

iT=ifum - — - — 

1RAII )= IRA ( 1ST ) 

I RA ( I ST » =1 T 

B8T*R9| 1ST ) _ , — - 

9R(IST»=FPm 

SP(t»=P?T 

R lf=R 1 1 1ST) ■ — 

ri n$T»=Pim 
91111 *R IT 

100 CONTINUE 

IP=t+l 

00 110 J=IP,N . ... 

CALL CMUI ( QTR,0T1 ,UNR ( J, I > ,UMI ( J, I) ,RP ( I) tRI I I.U — 

PR< J)=PR< JI-CTR 
P! ( Jl =RI ( J)~QT! 

110 CONTINUE 

120 CONTINUE 

CHECK. FOR INCONSISTENCY 

IFU0PTJ127, 127,125 

125 PRINT 4300 — 

PRINT 4C02, (PP(-J) » J =1 *N I 
4300 F0 ,V 1AT ( 5X, *R IGHT H4N0 SIDES/) 

4002 FORMAT (4C2O.10I i — - - 

127 DO 130 J*1,NEIG 
JJ=N- J ♦! 

CALL C AfiS( A CUN, Rfi ( JJ),Rl( J J 1 1. 

IFUOUM'EPR I 130,130,135 
130 CONTINUE 

NU“»N'JM*1 . . . - 

GO T 0 140 
135 1 1 P* I I ♦ l 

I E 1 0 ( L L , II P ) *NU» 

GO T 0 200 
140 CONTINUE 

- |p CONSISTENT THEN SACK SU8STITITE FOR GENERALIZED EIGENVECTOR 

GO T 0 60 . . - - 

2 00 CONTINUE 
RETURN 


265 

266 

267 

268 

269 

270 


RUNOUT INE C.MUL (A,q,Cl,C2,Dl,02) 

I"PLIC*.T Rt.'-L *3 (A 1-,0-YJ 

A=C 1*0 1- C 2*02 
8=CI*02+C2*C1 

return 

END 


271 


Vjo ; *! E C C I V t 



2 72 

273 

274 

275 

276 
277. 

278 

279 

280 
281 
28 2 

283 

284 

285 

286 
287 

2*e 

289 

290. 

291 

292 

293 

294 

295 
2S6 


257 

298 


299 

200 

3CI 

302 

203 

304 

205 

306 

307 
2C8 

209 

210 
311 
212 

213 

214 

315 

316 

2 17 
318 

215 
220 
321 

3 22 
223 

324 

325 
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IMPLICIT REAL *8(A-H,0-Y> 

6=01*01 +02 *02 
IF (6-1. 00-40 J 50. 100, 100 
50 PRINT 4000 

4000 FORMAT (5X, ‘0IVICE CHECK IN CDIV — DIVIDEND RETURNED* 

A = C l . 

8 = C 2 

00 TO 110 

_ 100 A=(C. 1*D1+C2*02)7E 

B= ( C2*01-C l*D2)/E 
110 CONTINUE 

RETURN _ .. 

END 

SUBROUTINE T~ AQS( A ,C1 ,C2) 

implicit real *8<a-h,o-y) 

A=0S0RT(C1#CI+C2*C2) 

RETURN .. 

END 

SUPROCTTNC SUM LENGTH, S.SH) __ 

IMPLICIT REAL *8< a'-H.J- vl 
D I MENS I CN WORK A ( 127), WORK AA( 128>,S< 12) 

EOU T VALENCE! JZ , ZJ ) 

EOL ! VALENCE ( WOP. K A ( 1) ,W0FKAM2) ) 

OfiL ZRO= 0 .OCO 

SU. M N=03l ZPC 

C 

C ZERO OUT THE ACCUMULATING ARRAY 

1000 DO 1010 L«t, 128 
1010 WCPKAA(L)=C8LZR0 

C DO JOHN'S ALGORITHM 

C 

DO 1050 1 = 1, LENGTH _ 

WORK = S ( I I 
LAST JZ=- l 

ZJ = WCRK . 

IF(ZJ) 1012, 1050, 1C13 

1012 ZJ=-ZJ 

1013 JZ=JZ/ 167772 16 

SUMN=WTRKA( JZ) tWORK 
WCSKA(JZ)=C8LZR0 

GO TO 1015 . 

1014 SUMN=S'IMN*w6rkA ( JZ ) 

WORKAI JZ ) = OBL ZRO 

1015 Z J = SUMN ... _ _ 

IFI ZJ >1016)1050, 1017 

1016 ZJ=ZJ 

1017 JZ = JZ / 1 677 72 16 _ 

IF( JZ-LASTJZ >1020,1030, 1020 

1020 LA ST JZ =J Z 

GO TO 1014 ... 

1030 WO s K A ( J Z ) = SUMN 
1050 CONTINUE 

SUM‘I=0BLZR0 . . 

DO 1060 L= It 128 
1060 SUM')= SUFN+WCRKA A (L) 

999 SM = S'JMN ... 

RETURN N 

END 



